      module bernstein
!
!     Module for manipulation of Bernstein polynomials.
!     Used in Bernstein-poly - a Bernstein polynomial collocation solution procedure for two-point BVP for ODE
!  
      use prec

      contains 

      function factorial(n) result(reslt)
!     finds factorial of n = n*n-1*...*2*1
!     It's not integer function to avoid overflow
      use lanczos_gamma
      implicit none

      real(dp) :: reslt
      integer, intent(in) :: n
      real(dp) :: t

      t = dble(n+1)
!     1.
!        reslt = lancz_gamma_g9(t)
!     2.
!        reslt = lancz_gamma_g7(t)
!     3.
        reslt = gamma(t) 

!     Recursive definition: keep (n < 15)  to avoid overflow in factorial calculation
!      if (n==1 .or.n==0) then
!        reslt = 1
!      else
!        reslt = real(n)*factorial(n-1)
!      endif

      end function factorial

      function binomial_coef(n,k)
!     finds binomial coefficient 'n over k'
      implicit none

      real(dp) :: binomial_coef     
      integer, intent(in) :: n,k

!     Using Gamma for factorial evaluation: n must be <= 170      
      binomial_coef = factorial(n) /( factorial(k)*factorial(n-k) )

!     Using property that ties binomial coefficient to Log_Gamma function: n must be <= 130
!     Log_Gamma available only in Fortran 2008 and later!!!     
!      binomial_coef = exp( log_gamma(dble(n+1)) - log_gamma(dble(k+1)) - log_gamma(dble(n-k+1)) )

      end function binomial_coef

!      recursive function bernstein_basis_fun_eval(i,n,a,b,t) result(reslt)
      function bernstein_basis_fun_eval(i,n,a,b,t) result(reslt)
!     Evaluates i-th Bernstein basis function of order n, at parameter value t, t \in [0,1]
      implicit none

      real(dp) :: reslt
      integer, intent(in) :: i,n
      real(dp), intent(in) :: a,b ! start and end of interval where funct. is defined
      real(dp), intent(in) :: t

!     Non-recursive definition:
      if (i==0 .and. n==0) then
        reslt = 1.0_dp
      elseif (i < 0 .or. i > n) then
        reslt = 0.0_dp ! important for derivative evaluation
      else
        reslt = binomial_coef(n,i) * (t-a)**dble(i) * (b-t)**dble(n-i) / (b-a)**dble(n)
       endif

!     Recursive definition:      
!      if (i==0 .and. n==0) then
!        reslt = 1.0_dp
!      elseif (i < 0 .or. i > n) then
!        reslt = 0.0_dp ! important for derivative evaluation
!      else
!        reslt = (b-t)/(b-a)*bernstein_basis_fun_eval(i,n-1,a,b,t) + t/(b-a)*bernstein_basis_fun_eval(i-1,n-1,a,b,t)
!      endif

      end function bernstein_basis_fun_eval

      function bernstein_basis_fun_maxval(i,n,a,b) result(reslt)
!     finds maximum of ith Bernstein basis function of order n 
      implicit none

      real(dp) :: reslt
      integer, intent(in) :: i,n
      real(dp), intent(in) :: a,b ! start and end of interval where funct. is defined
      real(dp) :: t

!     Maximum is : i**i * (1./n)**n * (n-i)**(n-i) * binomial_coef(n,i), and it occurs at:
      t = a + dble(i)/dble(n)*(b-a)

      reslt = bernstein_basis_fun_eval(i,n,a,b,t)

      end function bernstein_basis_fun_maxval

      function bernstein_basis_nodes(n,a,b) result(reslt)
!     Returns position of the maximum of i-th Bernstein basis function of order n,
!     stored in an array. It serves as the vector of nodes in [a,b] interval.
      implicit none

      integer, intent(in) :: n
      real(dp), dimension(0:n) :: reslt ! length (n+1)
      real(dp), intent(in) :: a,b ! start and end of interval where funct. is defined
      integer :: i

      do i = 0,n,1
      reslt(i) = a + (dble(i)/dble(n))*(b-a)
      enddo

      end function bernstein_basis_nodes

      function bernstein_basis_fun_mult(i,j,k,m,a,b,t) result(reslt)
!     Returns multiplication of two Bernstein bases Bij*Bkm
      implicit none

      real(dp) :: reslt
      integer, intent(in) :: i,j,k,m
      real(dp), intent(in) :: a,b ! start and end of interval where funct. is defined
      real(dp), intent(in) :: t

      reslt = (binomial_coef(i,j)*binomial_coef(m,k)/binomial_coef(j+m,i+k))*bernstein_basis_fun_eval(i+k,j+m,a,b,t)

      end function bernstein_basis_fun_mult

      recursive double precision function bernstein_basis_fun_derivative(k,i,n,a,b,t) result(reslt)
!     The derivative can be written as a combination of two polynomials of lower degree.
      implicit none

!      real(dp) :: reslt
      integer, intent(in) :: k    ! derivative order
      integer, intent(in) :: i,n
      real(dp), intent(in) :: a,b ! start and end of interval where funct. is defined
      real(dp), intent(in) :: t   ! parameter value a <= t <= b
      

      if (k > n) write(*,*) 'bernstein_basis_fun_derivative::Error: k > n!'

      if (k == 1) then
      reslt = (dble(n)/(b-a))*( bernstein_basis_fun_eval(i-1,n-1,a,b,t) - bernstein_basis_fun_eval(i,n-1,a,b,t) ) 
      else
      reslt = dble(n)/(b-a)*(bernstein_basis_fun_derivative(k-1,i-1,n-1,a,b,t)-bernstein_basis_fun_derivative(k-1,i,n-1,a,b,t))
      end if

      end function bernstein_basis_fun_derivative   

      function Bernstein_interpolant(beta,n,a,b,t) result(reslt)
!     Approximates function using Bernstein polynomials in the form f(x) = sum_i=0^i=n (beta(i)*B_i,n(x)); where beta(i)=f(i/n)
      implicit none

      real(dp) :: reslt
      real(dp), dimension(0:n), intent(in) :: beta ! Bernstein polynomial coefficients !0:n non-linear. n+1 linear
      integer, intent(in) :: n                     ! order of Bernstein polynomial
      real(dp), intent(in) :: a,b
      real(dp), intent(in) :: t ! position of evaluation
      integer :: i
      real(dp) :: suma

      suma = 0.0_dp
      do i=0,n,1 ! sum all the basis function at point x
        suma = suma + beta(i)*bernstein_basis_fun_eval(i,n,a,b,t)
      enddo
      reslt = suma

      end function Bernstein_interpolant


      end module bernstein
